Tile: Telepressure & Suicide in Secundary School Spain - Network Analysis: A Network Analysis Approach
Authors:
José-Manuel Guerra1,
Jesús Martínez-Calvo2,
Diego Diaz-Milanes2,3
Affiliation:
1 Department of Social Psychology, University of Seville, Seville, 41018, Spain
2 Department of Quantitative Methods, Universidad Loyola Andalucía, Sevilla, Spain
3 Health Research Institute, University of Canberra, Canberra, ACT, Australia
This R Markdown file contains the de-identified data and all R code used in the study, ensuring full transparency and reproducibility.
The workflow comprises ten steps:
Data import and preprocessing
Load the dataset, screen missingness/outliers, and prepare
variables.
Descriptive statistics
Summarize sociodemographics, inspect the distribution of variables, and
explore gender differences where relevant.
Undirected network modeling
Estimate a partial correlation (GGM) network of WTM
items using EBICglasso.
Network invariance testing
Compare network structure, global strength, and specific edges across
gender groups via the Network Comparison Test
(FDR-adjusted edgewise tests).
Directed network modeling
Learn a Bayesian Network (DAG) for potential
directional dependencies among WTM items with bootstrap
aggregation.
Predictive evaluation
Assess model predictability (e.g., \(R^2\), RMSE) via cross-validation and
compare performance across models.
In this section, the packages, function design, and data used are listed. Then, data processing—including variable selection and filtering by units of analysis—was performed.
# ---- Load libraries (already attached above; keep here if you knit this chunk standalone)
suppressPackageStartupMessages({
library(haven); library(tidyr); library(dplyr); library(corrplot); library(psych)
library(parameters); library(ggplot2); library(semTools)
library(tibble); library(purrr); library(effectsize); library(EGAnet); library(qgraph)
library(NetworkComparisonTest); library(bootnet); library(mgm); library(bnlearn); library(caret)
library(forcats); library(knitr); library(kableExtra)
})
# ---- Helper: safe version + loaded flag
safe_version <- function(pkg) {
v <- tryCatch(as.character(utils::packageVersion(pkg)), error = function(e) NA_character_)
if (is.na(v)) "" else v
}
is_loaded <- function(pkg) pkg %in% loadedNamespaces()
# ---- Packages + key functions actually used/reported
pkgs <- tibble::tibble(
Package = c(
"readxl","tidyr","dplyr","psych","semTools","ggplot2","effectsize",
"EGAnet","qgraph","NetworkComparisonTest","bootnet","mgm","bnlearn",
"caret","forcats","knitr","kableExtra","lavaan","semPlot","corrplot","haven","parameters"
),
Functions = c(
"`read_excel`",
"`pivot_longer`, `pivot_wider`, reshaping support",
"`select`, `mutate`, `filter`, `group_by`, `summarise`, `bind_rows`",
"`describe`, `polychoric`, `alpha`, `omega`, `splitHalf`",
"`mardiaSkew`, `mardiaKurtosis`, invariance helpers",
"`ggplot`, `aes`, `geom_line`, `labs`, `theme_bw`",
"`cohens_d`",
"`bootEGA`, `dimensionStability`",
"`qgraph`, `EBICglasso`",
"`NCT`, `plot`",
"`estimateNetwork`, `bootnet`, `corStability`",
"`mgm`, `predict`",
"`boot.strength`, `averaged.network`, `bn.fit`, `predict`",
"`confusionMatrix`",
"`fct_reorder`",
"`kable`",
"`kable_styling`, `add_header_above`",
"`cfa`, `fitMeasures`",
"`semPaths`",
"`corrplot`",
"`read_sav`, `as_factor`",
"`model_parameters`"
)
)
# ---- Add Version and Loaded columns
packages_info <- pkgs %>%
mutate(
Version = vapply(Package, safe_version, character(1)),
Loaded = ifelse(vapply(Package, is_loaded, logical(1)), "Yes", "No")
) %>%
select(Package, Version, Loaded, Functions)
# ---- Display table
knitr::kable(
packages_info,
caption = "Table 1. Packages used, versions, loaded status, and key functions"
) %>%
kableExtra::kable_styling(full_width = FALSE, position = "center")| Package | Version | Loaded | Functions |
|---|---|---|---|
| readxl | 1.4.3 | No |
read_excel
|
| tidyr | 1.3.1 | Yes |
pivot_longer, pivot_wider, reshaping support
|
| dplyr | 1.1.4 | Yes |
select, mutate, filter,
group_by, summarise, bind_rows
|
| psych | 2.4.1 | Yes |
describe, polychoric, alpha,
omega, splitHalf
|
| semTools | 0.5.6 | Yes |
mardiaSkew, mardiaKurtosis, invariance helpers
|
| ggplot2 | 3.5.2 | Yes |
ggplot, aes, geom_line,
labs, theme_bw
|
| effectsize | 0.8.9 | Yes |
cohens_d
|
| EGAnet | 2.0.4 | Yes |
bootEGA, dimensionStability
|
| qgraph | 1.9.8 | Yes |
qgraph, EBICglasso
|
| NetworkComparisonTest | 2.2.2 | Yes |
NCT, plot
|
| bootnet | 1.5.6 | Yes |
estimateNetwork, bootnet,
corStability
|
| mgm | 1.2.14 | Yes |
mgm, predict
|
| bnlearn | 4.9 | Yes |
boot.strength, averaged.network,
bn.fit, predict
|
| caret | 6.0.94 | Yes |
confusionMatrix
|
| forcats | 1.0.0 | Yes |
fct_reorder
|
| knitr | 1.45 | Yes |
kable
|
| kableExtra | 1.4.0 | Yes |
kable_styling, add_header_above
|
| lavaan | 0.6.17 | Yes |
cfa, fitMeasures
|
| semPlot | 1.1.6 | No |
semPaths
|
| corrplot | 0.92 | Yes |
corrplot
|
| haven | 2.5.4 | Yes |
read_sav, as_factor
|
| parameters | 0.22.1 | Yes |
model_parameters
|
Estamos utilizando 22 paquetes.
Design and generation of a formula for plotting centrality measures of two networks.
compareCentrality <- function(net1, net2,
include = c("Strength",
"Closeness",
"Betweenness",
"ExpectedInfluence",
"all",
"All"),
orderBy = c("Strength",
"Closeness",
"Betweenness",
"ExpectedInfluence"),
decreasing = T,
legendName = '',
net1Name = 'Network 1',
net2Name = 'Network 2'){
if(include == "All" | include == "all"){
include = c("Strength",
"Closeness",
"Betweenness",
"ExpectedInfluence")
}
df <- centralityTable(net1, net2) %>% filter(measure %in% include)
df %>%
mutate(graph = case_when(graph == 'graph 1' ~ net1Name,
graph == 'graph 2' ~ net2Name),
graph = as.factor(graph),
node = as.factor(node)) %>%
mutate(node = fct_reorder(node, value)) %>%
ggplot(aes(x = node, y = value, group = graph)) +
geom_line(aes(linetype = graph), size = 1) +
labs(x = '', y = '') +
scale_linetype_discrete(name = legendName) +
coord_flip() +
facet_grid(~measure) +
theme_bw()
}The dataset from the #### project (Ethical Committee: ####) was loaded, including only the minimum set of variables required for the present study.
Filters were applied to remove participants who did not complete the isntrument or indicate their gender, and to exclude cases with missing values:
tele_df <- dplyr::select(df_general,c(Sexo,EdadV1,V2,AUTO:TPB))
names(tele_df)[1:3] <- c("sex","age","academic_year")
tele_df <- as.data.frame(tele_df)
tele_df$sex <- as.factor(case_when(tele_df$sex == 1 ~ "Male",
tele_df$sex == 2 ~ "Female",
tele_df$sex == 3 ~ "I do not know"))
tele_df$academic_year <- as.factor(case_when(tele_df$academic_year == 1 ~ "First academic period",
tele_df$academic_year == 2 ~ "First academic period",
tele_df$academic_year == 3 ~ "Second academic period",
tele_df$academic_year == 4 ~ "Second academic period"))
tele_df$PS <- as.factor(case_when(tele_df$PS == 0 ~ "Passive-Agressive",
tele_df$PS == 1 ~ "Assertive",
tele_df$PS == 2 ~ "Agressive",
tele_df$PS == 3 ~ "Passive"))
tele_df$TPB <- as.factor(case_when(tele_df$TPB == 1 ~ "Low risk",
tele_df$TPB == 2 ~ "Medium risk",
tele_df$TPB == 3 ~ "High risk"))
tele_df <- dplyr::select(tele_df, c("sex","academic_year","PS","TPB","age","AUTO","HETERO",
"TP","SP","RS","CIUS","IDS","AB","VB","VCB","ACB","AF",
"AA","APS"))
initial_rows <- nrow(tele_df)
tele_df <- na.omit(tele_df)
print(paste("Initial number of rows:", initial_rows, "| Final number of rows:", nrow(tele_df)))## [1] "Initial number of rows: 820 | Final number of rows: 813"
Creation of dataframes tailored to meet the requirements of the forthcoming data analyses.
1. Dataframe for the analysis of individual items
The following table displays the first five rows of the set of items from the WTM instrument, used in the present study.
# Prepare the dataframe
network_variables <- dplyr::select(tele_df, age:APS)
pieColor <- rep("#EB9446", length(network_variables))
# Display first 5 rows with a formatted table
head(network_variables, 5) %>%
kable(caption = "Table 2. Set of items from the WTM") %>%
kable_styling(full_width = FALSE, position = "center")| age | AUTO | HETERO | TP | SP | RS | CIUS | IDS | AB | VB | VCB | ACB | AF | AA | APS |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 14 | 10 | 27 | 3.500000 | 2.75 | 3.333333 | 2.533333 | 2.4 | 1.333333 | 1.333333 | 1.000000 | 1.000000 | 2.25 | 3.00 | 2.75 |
| 13 | 55 | 7 | 2.500000 | 2.25 | 2.000000 | 2.866667 | 3.2 | 1.666667 | 2.333333 | 1.000000 | 2.000000 | 1.00 | 1.25 | 1.25 |
| 13 | 55 | 53 | 1.833333 | 3.25 | 1.666667 | 1.333333 | 1.2 | 1.000000 | 2.000000 | 1.666667 | 1.000000 | 4.00 | 1.50 | 4.00 |
| 14 | 50 | 53 | 2.333333 | 2.75 | 2.000000 | 1.800000 | 1.0 | 1.666667 | 2.000000 | 1.333333 | 1.333333 | 4.00 | 4.00 | 3.75 |
| 14 | 45 | 47 | 2.500000 | 2.00 | 1.666667 | 2.600000 | 3.8 | 1.666667 | 1.333333 | 1.000000 | 1.000000 | 1.25 | 2.75 | 2.75 |
2. Dataframe for exploring sex differences
The following table shows the first five rows of the dataset including the WTM item responses, the total score, and participants’ sex
# Select CEI items, total score, and gender
NV_sex <- dplyr::select(tele_df, c(sex,age:APS))
# Display first 5 rows with styling
head(NV_sex, 5) %>%
kable(caption = "Table 3. WTM items, total score, and gender variable") %>%
kable_styling(full_width = FALSE, position = "center")| sex | age | AUTO | HETERO | TP | SP | RS | CIUS | IDS | AB | VB | VCB | ACB | AF | AA | APS |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Female | 14 | 10 | 27 | 3.500000 | 2.75 | 3.333333 | 2.533333 | 2.4 | 1.333333 | 1.333333 | 1.000000 | 1.000000 | 2.25 | 3.00 | 2.75 |
| Female | 13 | 55 | 7 | 2.500000 | 2.25 | 2.000000 | 2.866667 | 3.2 | 1.666667 | 2.333333 | 1.000000 | 2.000000 | 1.00 | 1.25 | 1.25 |
| Female | 13 | 55 | 53 | 1.833333 | 3.25 | 1.666667 | 1.333333 | 1.2 | 1.000000 | 2.000000 | 1.666667 | 1.000000 | 4.00 | 1.50 | 4.00 |
| Male | 14 | 50 | 53 | 2.333333 | 2.75 | 2.000000 | 1.800000 | 1.0 | 1.666667 | 2.000000 | 1.333333 | 1.333333 | 4.00 | 4.00 | 3.75 |
| Female | 14 | 45 | 47 | 2.500000 | 2.00 | 1.666667 | 2.600000 | 3.8 | 1.666667 | 1.333333 | 1.000000 | 1.000000 | 1.25 | 2.75 | 2.75 |
2. Dataframe for exploring academic year differences
The following table shows the first five rows of the dataset including the WTM item responses, the total score, and participants’ academic year
NV_academic_year <- dplyr::select(tele_df, c(academic_year,age:APS))
# Display first 5 rows with styling
head(NV_academic_year, 5) %>%
kable(caption = "Table 3. WTM items, total score, and gender variable") %>%
kable_styling(full_width = FALSE, position = "center")| academic_year | age | AUTO | HETERO | TP | SP | RS | CIUS | IDS | AB | VB | VCB | ACB | AF | AA | APS |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| First academic period | 14 | 10 | 27 | 3.500000 | 2.75 | 3.333333 | 2.533333 | 2.4 | 1.333333 | 1.333333 | 1.000000 | 1.000000 | 2.25 | 3.00 | 2.75 |
| First academic period | 13 | 55 | 7 | 2.500000 | 2.25 | 2.000000 | 2.866667 | 3.2 | 1.666667 | 2.333333 | 1.000000 | 2.000000 | 1.00 | 1.25 | 1.25 |
| First academic period | 13 | 55 | 53 | 1.833333 | 3.25 | 1.666667 | 1.333333 | 1.2 | 1.000000 | 2.000000 | 1.666667 | 1.000000 | 4.00 | 1.50 | 4.00 |
| First academic period | 14 | 50 | 53 | 2.333333 | 2.75 | 2.000000 | 1.800000 | 1.0 | 1.666667 | 2.000000 | 1.333333 | 1.333333 | 4.00 | 4.00 | 3.75 |
| First academic period | 14 | 45 | 47 | 2.500000 | 2.00 | 1.666667 | 2.600000 | 3.8 | 1.666667 | 1.333333 | 1.000000 | 1.000000 | 1.25 | 2.75 | 2.75 |
Participants distribution by sex and age:
# sex distribution
sex_table <- table(tele_df$sex)
sex_percent <- round(prop.table(sex_table) * 100, 2)
sex_df <- data.frame(
sex = names(sex_table),
Frequency = round(as.numeric(sex_table), 2),
Percentage = paste0(sex_percent, "%")
)
# Age summary
df_age <- tele_df[!is.na(tele_df$age), ]
age_summary <- summary(df_age$age)
age_sd <- round(sd(df_age$age), 2)
age_df <- data.frame(
Statistic = names(age_summary),
Value = round(as.numeric(age_summary), 2)
) %>%
bind_rows(data.frame(Statistic = "Standard Deviation", Value = age_sd))
# Display tables
kable(sex_df, caption = "Table 5. sex distribution (frequency and %)") %>%
kable_styling(full_width = FALSE, position = "center")| sex | Frequency | Percentage |
|---|---|---|
| Female | 400 | 49.2% |
| I do not know | 10 | 1.23% |
| Male | 403 | 49.57% |
kable(age_df, caption = "Table 6. Age summary statistics") %>%
kable_styling(full_width = FALSE, position = "center")| Statistic | Value |
|---|---|
| Min. | 12.00 |
| 1st Qu. | 13.00 |
| Median | 14.00 |
| Mean | 13.84 |
| 3rd Qu. | 15.00 |
| Max. | 18.00 |
| Standard Deviation | 1.24 |
Univariate descriptive analysis of the virables included in the analysis:
# Descriptive statistics
descrgen <- describe(network_variables)
# Build dataframe with selected stats
lgen <- list(c(1:15), descrgen$mean, descrgen$sd, descrgen$min,
descrgen$max, descrgen$skew, descrgen$kurtosis)
lgen <- as.data.frame(lgen)
floor_pct <- sapply(network_variables, function(x) mean(x == min(x, na.rm = TRUE), na.rm = TRUE) * 100)
ceiling_pct <- sapply(network_variables, function(x) mean(x == max(x, na.rm = TRUE), na.rm = TRUE) * 100)
lgen <- cbind(lgen,floor_pct)
lgen <- cbind(lgen,ceiling_pct)
# Round and rename columns
lgen <- lgen %>%
mutate_if(is.numeric, round, digits = 3)
colnames(lgen) <- c("Item", "Mean", "SD", "Min", "Max", "Skew", "Kurtosis", "Floor%", "Ceiling%")
# Display table
kable(lgen, caption = "Table 7. Descriptive statistics for WTM items") %>%
kable_styling(full_width = FALSE, position = "center")| Item | Mean | SD | Min | Max | Skew | Kurtosis | Floor% | Ceiling% | |
|---|---|---|---|---|---|---|---|---|---|
| age | 1 | 13.835 | 1.240 | 12 | 18.000 | 0.303 | -0.415 | 15.744 | 0.246 |
| AUTO | 2 | 67.583 | 24.830 | 0 | 100.000 | -0.574 | -0.454 | 0.861 | 11.808 |
| HETERO | 3 | 59.513 | 26.603 | 0 | 100.000 | -0.148 | -0.871 | 1.353 | 11.070 |
| TP | 4 | 2.486 | 0.822 | 1 | 5.000 | 0.298 | -0.229 | 4.182 | 0.246 |
| SP | 5 | 2.230 | 0.804 | 1 | 4.000 | 0.237 | -0.770 | 10.578 | 2.583 |
| RS | 6 | 2.211 | 0.893 | 1 | 4.000 | 0.290 | -0.984 | 15.867 | 4.797 |
| CIUS | 7 | 2.427 | 0.680 | 1 | 4.571 | 0.324 | -0.101 | 0.984 | 0.123 |
| IDS | 8 | 1.490 | 0.825 | 1 | 5.000 | 2.105 | 4.044 | 52.399 | 0.492 |
| AB | 9 | 1.484 | 0.700 | 1 | 5.000 | 2.257 | 6.099 | 44.895 | 0.615 |
| VB | 10 | 1.860 | 0.899 | 1 | 5.000 | 1.362 | 1.614 | 24.477 | 1.353 |
| VCB | 11 | 1.394 | 0.705 | 1 | 5.000 | 2.642 | 8.209 | 59.410 | 0.984 |
| ACB | 12 | 1.280 | 0.595 | 1 | 5.000 | 3.380 | 13.940 | 67.651 | 0.492 |
| AF | 13 | 3.354 | 0.765 | 1 | 4.000 | -1.316 | 1.046 | 2.091 | 36.039 |
| AA | 14 | 3.230 | 0.803 | 1 | 4.000 | -0.762 | -0.554 | 0.861 | 34.317 |
| APS | 15 | 3.290 | 0.745 | 1 | 4.000 | -0.926 | 0.056 | 0.984 | 33.948 |
Assessment of multivariate normality assumptions based on Mardia’s tests for skewness and kurtosis, applied to the set of items included in the WTM instrument:
# Run Mardia's tests
skew_result <- mardiaSkew(network_variables)
kurt_result <- mardiaKurtosis(network_variables)
# Round and extract values
mardia_df <- data.frame(
Test = c("Mardia Skewness", "Mardia Kurtosis"),
Statistic = round(c(skew_result[2], kurt_result[2]), 3),
p_value = round(c(skew_result[4], kurt_result[3]), 3)
)
row.names(mardia_df) <- 1:2
# Show table
kable(mardia_df, caption = "Table 8: Mardia's Tests for Multivariate Normality")| Test | Statistic | p_value |
|---|---|---|
| Mardia Skewness | 5494.117 | 0 |
| Mardia Kurtosis | 48.142 | 0 |
# Lista de ítems y etiquetas
items <- names(network_variables) #PARA MODIFICAR
item_labels <- c(
"Age",
"Autoassertivity",
"Heteroassertivity",
"Telepressure",
"Supervision",
"Restriction",
"Problematic use of internet",
"Suicide ideation",
"Agression bullying",
"Victim bullying",
"Victim ciberbullying",
"Agressor ciberbullying",
"Family support",
"Friends support",
"Significant Person support"
)
NV_sex <- NV_sex[NV_sex$sex != "I do not know",]
# Crear tabla de resultados
results <- lapply(seq_along(items), function(i) {
item <- items[i]
label <- item_labels[i]
# Obtener estadísticas por grupo
stats <- NV_sex %>%
group_by(sex) %>%
summarise(
mean = round(mean(.data[[item]], na.rm = TRUE), 3),
sd = round(sd(.data[[item]], na.rm = TRUE), 3)
)
male_val <- paste0(stats$mean[stats$sex == "Male"], " (", stats$sd[stats$sex == "Male"], ")")
female_val <- paste0(stats$mean[stats$sex == "Female"], " (", stats$sd[stats$sex == "Female"], ")")
# t-test
t_result <- t.test(as.formula(paste(item, "~ sex")), data = NV_sex)
t_stat <- round(t_result$statistic, 3)
df <- round(t_result$parameter, 3)
p_val <- round(t_result$p.value,3)
# Efecto (Cohen's d)
d_val <- cohens_d(as.formula(paste(item, "~ sex")), data = NV_sex)$Cohens_d
d_val <- round(d_val, 3)
# Combinar resultados
data.frame(
Item = label,
Male = male_val,
Female = female_val,
T_statistic = paste0(t_stat, " (", df, ")"),
p_value = p_val,
Effect_size = d_val
)
})
# Unir todos los ítems
final_table <- bind_rows(results)
# Mostrar tabla
kable(final_table, caption = "Table 9. Mean (SD) by sex, t-test Results, p-value and Effect Sizes for numeric variables")| Item | Male | Female | T_statistic | p_value | Effect_size |
|---|---|---|---|---|---|
| Age | 13.878 (1.251) | 13.773 (1.214) | -1.217 (800.604) | 0.224 | -0.086 |
| Autoassertivity | 71.998 (23.962) | 63.163 (24.869) | -5.126 (799.408) | 0.000 | -0.362 |
| Heteroassertivity | 65.089 (26.07) | 53.8 (25.832) | -6.164 (800.998) | 0.000 | -0.435 |
| Telepressure | 2.41 (0.827) | 2.564 (0.807) | 2.667 (800.759) | 0.008 | 0.188 |
| Supervision | 2.181 (0.799) | 2.294 (0.802) | 2.008 (800.885) | 0.045 | 0.142 |
| Restriction | 2.201 (0.91) | 2.232 (0.877) | 0.493 (800.316) | 0.622 | 0.035 |
| Problematic use of internet | 2.364 (0.638) | 2.48 (0.708) | 2.437 (791.239) | 0.015 | 0.172 |
| Suicide ideation | 1.363 (0.747) | 1.597 (0.856) | 4.138 (785.192) | 0.000 | 0.292 |
| Agression bullying | 1.589 (0.747) | 1.355 (0.581) | -4.964 (757.796) | 0.000 | -0.350 |
| Victim bullying | 1.946 (0.959) | 1.76 (0.804) | -2.975 (779.264) | 0.003 | -0.210 |
| Victim ciberbullying | 1.403 (0.727) | 1.367 (0.644) | -0.729 (791.029) | 0.466 | -0.051 |
| Agressor ciberbullying | 1.319 (0.666) | 1.228 (0.459) | -2.264 (714.003) | 0.024 | -0.160 |
| Family support | 3.459 (0.679) | 3.268 (0.819) | -3.611 (772.53) | 0.000 | -0.255 |
| Friends support | 3.267 (0.776) | 3.21 (0.812) | -1.012 (798.818) | 0.312 | -0.071 |
| Significant Person support | 3.316 (0.727) | 3.281 (0.74) | -0.666 (800.519) | 0.505 | -0.047 |
# Lista de ítems y etiquetas
items <- names(network_variables) #PARA MODIFICAR
item_labels <- c(
"Age",
"Autoassertivity",
"Heteroassertivity",
"Telepressure",
"Supervision",
"Restriction",
"Problematic use of internet",
"Suicide ideation",
"Agression bullying",
"Victim bullying",
"Victim ciberbullying",
"Agressor ciberbullying",
"Family support",
"Friends support",
"Significant Person support"
)
# Crear tabla de resultados
results <- lapply(seq_along(items), function(i) {
item <- items[i]
label <- item_labels[i]
# Obtener estadísticas por grupo
stats <- NV_academic_year %>%
group_by(academic_year) %>%
summarise(
mean = round(mean(.data[[item]], na.rm = TRUE), 3),
sd = round(sd(.data[[item]], na.rm = TRUE), 3)
)
male_val <- paste0(stats$mean[stats$academic_year == "First academic period"], " (", stats$sd[stats$academic_year == "First academic period"], ")")
female_val <- paste0(stats$mean[stats$academic_year == "Second academic period"], " (", stats$sd[stats$academic_year == "Second academic period"], ")")
# t-test
t_result <- t.test(as.formula(paste(item, "~ academic_year")), data = NV_academic_year)
t_stat <- round(t_result$statistic, 3)
df <- round(t_result$parameter, 3)
p_val <- round(t_result$p.value,3)
# Efecto (Cohen's d)
d_val <- cohens_d(as.formula(paste(item, "~ academic_year")), data = NV_academic_year)$Cohens_d
d_val <- round(d_val, 3)
# Combinar resultados
data.frame(
Item = label,
First_Period = male_val,
Second_Period = female_val,
T_statistic = paste0(t_stat, " (", df, ")"),
p_value = p_val,
Effect_size = d_val
)
})
# Unir todos los ítems
final_table <- bind_rows(results)
# Mostrar tabla
kable(final_table, caption = "Table 10. Mean (SD) by academic year, t-test Results, p-value and Effect Sizes for numeric variables")| Item | First_Period | Second_Period | T_statistic | p_value | Effect_size |
|---|---|---|---|---|---|
| Age | 13.03 (0.817) | 14.911 (0.822) | -32.377 (745.26) | 0.000 | -2.297 |
| Autoassertivity | 68.548 (24.376) | 66.293 (25.403) | 1.274 (730.568) | 0.203 | 0.091 |
| Heteroassertivity | 60.686 (26.028) | 57.945 (27.312) | 1.444 (727.618) | 0.149 | 0.103 |
| Telepressure | 2.434 (0.813) | 2.556 (0.83) | -2.101 (739.134) | 0.036 | -0.149 |
| Supervision | 2.439 (0.795) | 1.949 (0.729) | 9.12 (778.592) | 0.000 | 0.638 |
| Restriction | 2.42 (0.874) | 1.932 (0.842) | 8.037 (761.665) | 0.000 | 0.567 |
| Problematic use of internet | 2.392 (0.669) | 2.472 (0.692) | -1.66 (733.82) | 0.097 | -0.118 |
| Suicide ideation | 1.526 (0.854) | 1.441 (0.782) | 1.478 (778.999) | 0.140 | 0.103 |
| Agression bullying | 1.453 (0.661) | 1.526 (0.748) | -1.435 (693.834) | 0.152 | -0.104 |
| Victim bullying | 1.883 (0.893) | 1.83 (0.908) | 0.834 (740.582) | 0.405 | 0.059 |
| Victim ciberbullying | 1.351 (0.611) | 1.452 (0.811) | -1.957 (620.654) | 0.051 | -0.144 |
| Agressor ciberbullying | 1.212 (0.477) | 1.371 (0.715) | -3.591 (569.714) | 0.000 | -0.269 |
| Family support | 3.375 (0.761) | 3.326 (0.77) | 0.907 (742.657) | 0.364 | 0.064 |
| Friends support | 3.169 (0.805) | 3.31 (0.793) | -2.492 (753.129) | 0.013 | -0.176 |
| Significant Person support | 3.282 (0.731) | 3.301 (0.763) | -0.363 (729.875) | 0.717 | -0.026 |
A bootstrap exploratory graph analysis (EGA) was performed to estimate the number, structure, and stability of variable communities (i.e., network dimensions) based on their structural consistency. In this context, a community refers to a cluster of nodes that exhibit strong interconnections through edges.
The analysis facilitates the identification and delineation of such communities, as well as the allocation of individual items to their corresponding dimensions.
set.seed(123)
communities <- bootEGA(network_variables,
model = "glasso",
type = "resampling",
iter = 500,
typicalStructure = TRUE,
plot.typicalStructure = TRUE)Figure 1. Community Estimation
Items were retained within their assigned communities only when their stability coefficients exceeded 0.70. Items falling below this threshold were excluded due to their potential to undermine the structural integrity of the corresponding dimension.
Figure 2. Community Stability
# Definir los grupos (según índices de df_net)
group_list <- list(
"Parenting Control" = c(1,5,6),
"Social-Digital Use" = c(2,3,4,7),
"Bullying" = c(9:12),
"Social Support" = c(13:15),
"Suicide"= 8
)
TP_glasso <- qgraph::qgraph(
cor_auto(network_variables),
graph = "glasso",
layout = "spring",
groups = group_list,
#color = group_colors,
palette = "colorblind",
nodeNames = item_labels,
labels = TRUE,
threshold = TRUE,
sampleSize = nrow(network_variables),
vsize = 5,
esize = 5,
asize = 6,
border.width = 1,
border.color = "black",
unCol = "black",
theme = "colorblind",
legend = TRUE, # <‑‑ leyenda automática
legend.cex = 0.8, # tamaño texto de la leyenda
layoutScale = c(0.8, 0.8), # encoge ligeramente la red
rescale = TRUE
)Figure 3. Gaussian Graphical Model Estimation
possible_edges <- (ncol(network_variables)*(ncol(network_variables)-1))/2
real_edges <- length(TP_glasso$Edgelist$weight)
density <- round(real_edges/possible_edges*100, digits = 2)
# Print explanation and result
cat("Number of potential edges:", possible_edges, "\n")## Number of potential edges: 105
## Number of real edges: 19
cat("Network density (percentage of actual edges out of all possible undirected edges):", density, "%\n")## Network density (percentage of actual edges out of all possible undirected edges): 18.1 %
# Extract edge weights
edge_weights <- TP_glasso$Edgelist$weight
# Calculate summary statistics
mean_weight <- round(mean(edge_weights), digits = 3)
sd_weight <- round(sd(edge_weights), digits = 3)
max_weight <- round(max(edge_weights), digits = 3)
min_weight <- round(min(edge_weights), digits = 3)
# Print explanation and results clearly
cat("Standard Edge Weight \n")## Standard Edge Weight
## The average edge weight was: 0.178 (SD: 0.298 )
## Minimum weight: -0.245 | Maximum weight: 0.574
# Calculate summary statistics
mean_weight <- round(mean(abs(edge_weights)), digits = 3)
sd_weight <- round(sd(abs(edge_weights)), digits = 3)
max_weight <- round(max(abs(edge_weights)), digits = 3)
min_weight <- round(min(abs(edge_weights)), digits = 3)
# Print explanation and results clearly
cat("\n")## Absolute Values
## The average edge weight was: 0.298 (SD: 0.169 )
## Minimum weight: 0.105 | Maximum weight: 0.574
TP_matrix <- as.matrix(network_variables)
invisible(capture.output({SV_mgm <-mgm(
data = TP_matrix,
type = rep("g", length(network_variables)),
level = rep(1, length(network_variables)),
k = 2,
verbatim = TRUE,
warnings = TRUE
)
}))
#Predice la varianza explicada de cada nodo y otros parametros
SV_predic <- predict(SV_mgm, network_variables, error.continuous = "VarExpl")
TP_glasso_predicted <- qgraph::qgraph(
cor_auto(network_variables),
graph = "glasso",
layout = "spring",
groups = group_list,
#color = group_colors,
palette = "colorblind",
nodeNames = item_labels,
labels = TRUE,
threshold = TRUE,
sampleSize = nrow(network_variables),
vsize = 5,
esize = 5,
asize = 6,
border.width = 1,
border.color = "black",
unCol = "black",
theme = "colorblind",
pie = SV_predic$errors$R2,
pieColor = pieColor,
legend = TRUE, # <‑‑ leyenda automática
legend.cex = 0.8, # tamaño texto de la leyenda
layoutScale = c(0.8, 0.8), # encoge ligeramente la red
rescale = TRUE
)Figure 4. Gaussian Graphical Model Characterization
centrality_TP <- centralityPlot(TP_glasso, include =
c("Betweenness","Closeness","Strength","ExpectedInfluence"),
orderBy ="Betweenness", scale = "z-scores")## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
Figure 5. Centrality Measures Estimation
# 1. Estimar la red con EBICglasso
glasso_net <- estimateNetwork(network_variables,
default = "EBICglasso",
corMethod = "cor_auto",
verbose = FALSE)
# 2. Bootstrap de estabilidad (case-dropping)
set.seed(2025)
boot_results <- invisible(bootnet(
estimateNetwork(network_variables, default = "EBICglasso", corMethod = "cor_auto"),
type = "case",
nBoots = 200,
nCores = 1, # importante para ver mensajes
statistics= c("strength", "closeness", "betweenness","expectedInfluence"), # evita centralidades inestables
verbose = FALSE,
maxErrors = 1000
))
# 3. Plot de estabilidad de centralidad (bootstrap case-dropping)
plot(boot_results,
statistics = c("strength", "closeness", "betweenness","expectedInfluence"),
labels = TRUE, order = "sample", legend = TRUE,
color = "darkblue", line = TRUE,
main = "Bootstrap Centrality Stability (Case Dropping)")Figure 6. Centrality Measures Stability
# 4. Coeficientes de estabilidad de centralidad
cs_coeffs <- corStability(boot_results, verbose = FALSE)
# Crear tabla
cs_table <- data.frame(
Centrality = c("Betweenness", "Closeness", "Strength", "ExpectedInfluence"),
CS_Coefficient = round(c(cs_coeffs[1], cs_coeffs[2], cs_coeffs[4], cs_coeffs[3]), 3)
)
row.names(cs_table) <- 1:4
# Mostrar como tabla en R Markdown
knitr::kable(cs_table, caption = "Table 10: Centrality Stability Coefficients (CS)",
align = c("r","r","r","r")) %>%
kableExtra::kable_styling(full_width = FALSE, position = "center", bootstrap_options = c("striped","hover"))| Centrality | CS_Coefficient |
|---|---|
| Betweenness | 0.594 |
| Closeness | 0.594 |
| Strength | 0.750 |
| ExpectedInfluence | 0.750 |
TP_male <- NV_sex[NV_sex$sex == "Male",-1]
TP_female <- NV_sex[NV_sex$sex == "Female",-1]
set.seed(1234)
EN_male = estimateNetwork(TP_male, default = "EBICglasso", corMethod = "cor_auto")
EN_female = estimateNetwork(TP_female, default = "EBICglasso", corMethod = "cor_auto")set.seed(1234)
invisible(capture.output({
res <- NCT(
TP_female,
TP_male,
p.adjust.methods = "BH",
binary.data = FALSE,
test.edges = TRUE,
edges = "all",
it = 1000
)
}))
female_str <- round(res$glstrinv.sep[1], 3)
male_str <- round(res$glstrinv.sep[2], 3)
p_str <- ifelse(res$glstrinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$glstrinv.pval, 3)))
dif_net <- round(res$nwinv.real, 3)
p_net <- ifelse(res$nwinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$nwinv.pval, 3)))
p_edges <- paste(
dplyr::filter(res$einv.pvals, `p-value` < 0.05) %>%
dplyr::mutate(pair = paste(Var1, "-", Var2)) %>%
dplyr::pull(pair),
collapse = ", "
)
pvalue_edges <- paste(
dplyr::filter(res$einv.pvals, `p-value` < 0.05)
)
pvalue_I4_I6 <- round(as.numeric(pvalue_edges[3]),3)
statistic_I4_I6 <- round(as.numeric(pvalue_edges[4]),3)In terms of global strength, the female network showed a similary overall connectivity (S = 5.893) compared to the male network (S = 5.745). See Figure 7.
fmt_num <- function(x, d = 3) {
if (is.numeric(x) && length(x) == 1 && is.finite(x)) sprintf(paste0("%.", d, "f"), x) else "NA"
}
fmt_p <- function(p, d = 3, eps = 1e-3) {
if (is.numeric(p) && length(p) == 1 && is.finite(p)) format.pval(p, digits = d, eps = eps) else "NA"
}
plot(res, what="strength")Figure 7. P-values for differences in network strength between sexs
However, the network structure invariance test indicated a difference of 0.242 between the female and male networks. This difference was statistically significant (p = 0.049), indicating that the overall network structures differ (Figure 8). Follow-up edgewise tests, with p-values adjusted for multiple comparisons using the Benjamini–Hochberg procedure, showed that the behavior of the edge Item4–Item6 was statistically significant (E = NA, p = NA).
Figure 8. P-values for differences in network structure between sexs
## [1] ""
The Network Comparison Test did not identify any statistically significant differences in edge strength between the UN/NW and OWO networks after applying the Benjamini-Hochberg procedure to correct p-values for multiple comparisons. The significance level was set at the conventional threshold of 0.05.
Female’s Network
# --- FEMALE GROUP ---
# Convert data to matrix
TP_female_matrix <- as.matrix(TP_female)
# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
TP_mgm_female <- mgm(
data = TP_female_matrix,
type = rep("g", length(TP_female)),
level = rep(1, length(TP_female)),
k = 2,
warnings = TRUE,
verbatim = TRUE
)
}))
# Predict explained variance (R²) and other node-level metrics
TP_predic_female <- predict(
TP_mgm_female,
TP_female,
error.continuous = "VarExpl"
)
TP_glasso_predicted_female <- qgraph::qgraph(
cor_auto(TP_female),
graph = "glasso",
layout = TP_glasso_predicted$layout,
groups = group_list,
#color = group_colors,
palette = "colorblind",
nodeNames = item_labels,
labels = TRUE,
threshold = TRUE,
sampleSize = nrow(network_variables),
vsize = 5,
esize = 5,
asize = 6,
border.width = 1,
border.color = "black",
unCol = "black",
theme = "colorblind",
pie = SV_predic$errors$R2,
pieColor = pieColor,
legend = TRUE, # <‑‑ leyenda automática
legend.cex = 0.8, # tamaño texto de la leyenda
layoutScale = c(0.8, 0.8), # encoge ligeramente la red
rescale = TRUE
)Figure 9. Networks of Female Participants
Male’s network
# --- MALE GROUP ---
# Convert data to matrix
TP_male_matrix <- as.matrix(TP_male)
# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
TP_mgm_male <- mgm(
data = TP_male_matrix,
type = rep("g", length(TP_male)),
level = rep(1, length(TP_male)),
k = 2,
warnings = TRUE,
verbatim = TRUE
)
}))
# Predict explained variance (R²) and other node-level metrics
TP_predic_male <- predict(
TP_mgm_male,
TP_male,
error.continuous = "VarExpl"
)
# Plot the network with R² represented as pie charts on nodes
TP_glasso_predicted_male <- qgraph::qgraph(
cor_auto(TP_male),
graph = "glasso",
layout = TP_glasso_predicted$layout,
groups = group_list,
#color = group_colors,
palette = "colorblind",
nodeNames = item_labels,
labels = TRUE,
threshold = TRUE,
sampleSize = nrow(network_variables),
vsize = 5,
esize = 5,
asize = 6,
border.width = 1,
border.color = "black",
unCol = "black",
theme = "colorblind",
pie = SV_predic$errors$R2,
pieColor = pieColor,
legend = TRUE, # <‑‑ leyenda automática
legend.cex = 0.8, # tamaño texto de la leyenda
layoutScale = c(0.8, 0.8), # encoge ligeramente la red
rescale = TRUE
)Figure 10. Networks of Male Participants
To evaluate the invariance of the centrality measures network across sexs.
# Compare centralities between networks
network_comparison <- compareCentrality(
TP_glasso_predicted_female,
TP_glasso_predicted_male,
include = "all",
legendName = "Centrality Measures by sex",
net1Name = "Female",
net2Name = "Male"
)
plot(network_comparison)Figure 11. Centrality comparison between female and male networks
set.seed(1234)
invisible(capture.output({
nct_test_centrality <- NCT(
EN_female,
EN_male,
it = 1000,
test.centrality = TRUE,
p.adjust.methods = "BH",
centrality = c("closeness", "betweenness", "strength", "expectedInfluence"),
)
}))
# Extract p-values and convert to data frame
centrality_pvals <- as.data.frame(nct_test_centrality$diffcen.pval)
# Round to 3 decimals (and adjust "< 0.001" if you want a more elegant output)
centrality_pvals_rounded <- round(centrality_pvals, 3)
# Optional: add row names as a column
centrality_pvals_rounded$Node <- rownames(centrality_pvals_rounded)
centrality_pvals_rounded <- centrality_pvals_rounded[, c("Node", setdiff(names(centrality_pvals_rounded), "Node"))]
row.names(centrality_pvals_rounded) <- 1:15
# Display as table
kable(centrality_pvals_rounded, caption = "Table 11. P-values for Centrality Differences (NCT Test)") %>%
kable_styling(full_width = FALSE, position = "center")| Node | closeness | betweenness | strength | expectedInfluence |
|---|---|---|---|---|
| age | 0.673 | 1.000 | 0.885 | 0.885 |
| AUTO | 0.587 | 1.000 | 1.000 | 1.000 |
| HETERO | 0.587 | 1.000 | 0.742 | 1.000 |
| TP | 0.587 | 0.885 | 1.000 | 0.587 |
| SP | 0.687 | 1.000 | 0.587 | 0.587 |
| RS | 0.687 | 1.000 | 0.673 | 1.000 |
| CIUS | 0.587 | 0.587 | 0.060 | 0.885 |
| IDS | 0.687 | 1.000 | 0.673 | 0.673 |
| AB | 0.603 | 0.901 | 1.000 | 0.673 |
| VB | 0.741 | 1.000 | 0.776 | 0.587 |
| VCB | 0.741 | 1.000 | 0.645 | 0.587 |
| ACB | 0.687 | 1.000 | 0.587 | 0.741 |
| AF | 0.685 | 1.000 | 0.587 | 0.587 |
| AA | 0.776 | 1.000 | 0.687 | 1.000 |
| APS | 0.885 | 1.000 | 0.587 | 1.000 |
TP_first <- NV_academic_year[NV_academic_year$academic_year == "First academic period",-1]
TP_second <- NV_academic_year[NV_academic_year$academic_year == "Second academic period",-1]
set.seed(12)
EN_first = estimateNetwork(TP_first, default = "EBICglasso", corMethod = "cor_auto")
EN_second = estimateNetwork(TP_second, default = "EBICglasso", corMethod = "cor_auto")set.seed(1234)
invisible(capture.output({
res <- NCT(
TP_second,
TP_first,
p.adjust.methods = "BH",
binary.data = FALSE,
test.edges = TRUE,
edges = "all",
it = 1000
)
}))
second_str <- round(res$glstrinv.sep[1], 3)
first_str <- round(res$glstrinv.sep[2], 3)
p_str <- ifelse(res$glstrinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$glstrinv.pval, 3)))
dif_net <- round(res$nwinv.real, 3)
p_net <- ifelse(res$nwinv.pval < 0.001, "p < 0.001", paste0("p = ", round(res$nwinv.pval, 3)))
p_edges <- paste(
dplyr::filter(res$einv.pvals, `p-value` < 0.05) %>%
dplyr::mutate(pair = paste(Var1, "-", Var2)) %>%
dplyr::pull(pair),
collapse = ", "
)
pvalue_edges <- paste(
dplyr::filter(res$einv.pvals, `p-value` < 0.05)
)
pvalue_I4_I6 <- round(as.numeric(pvalue_edges[3]),3)
statistic_I4_I6 <- round(as.numeric(pvalue_edges[4]),3)In terms of global strength, the second network showed a similary overall connectivity (S = 5.243) compared to the first network (S = 5.779). See Figure 7.
fmt_num <- function(x, d = 3) {
if (is.numeric(x) && length(x) == 1 && is.finite(x)) sprintf(paste0("%.", d, "f"), x) else "NA"
}
fmt_p <- function(p, d = 3, eps = 1e-3) {
if (is.numeric(p) && length(p) == 1 && is.finite(p)) format.pval(p, digits = d, eps = eps) else "NA"
}
plot(res, what="strength")Figure 7. P-values for differences in network strength between sexs
However, the network structure invariance test indicated a difference of 0.146 between the second and first networks. This difference was statistically significant (p = 0.512), indicating that the overall network structures differ (Figure 8). Follow-up edgewise tests, with p-values adjusted for multiple comparisons using the Benjamini–Hochberg procedure, showed that the behavior of the edge Item4–Item6 was statistically significant (E = NA, p = NA).
Figure 8. P-values for differences in network structure between sexs
## [1] ""
The Network Comparison Test did not identify any statistically significant differences in edge strength between the UN/NW and OWO networks after applying the Benjamini-Hochberg procedure to correct p-values for multiple comparisons. The significance level was set at the conventional threshold of 0.05.
second’s Network
# --- second GROUP ---
# Convert data to matrix
TP_second_matrix <- as.matrix(TP_second)
# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
TP_mgm_second <- mgm(
data = TP_second_matrix,
type = rep("g", length(TP_second)),
level = rep(1, length(TP_second)),
k = 2,
warnings = TRUE,
verbatim = TRUE
)
}))
# Predict explained variance (R²) and other node-level metrics
TP_predic_second <- predict(
TP_mgm_second,
TP_second,
error.continuous = "VarExpl"
)
TP_glasso_predicted_second <- qgraph::qgraph(
cor_auto(TP_second),
graph = "glasso",
layout = TP_glasso_predicted$layout,
groups = group_list,
#color = group_colors,
palette = "colorblind",
nodeNames = item_labels,
labels = TRUE,
threshold = TRUE,
sampleSize = nrow(network_variables),
vsize = 5,
esize = 5,
asize = 6,
border.width = 1,
border.color = "black",
unCol = "black",
theme = "colorblind",
pie = SV_predic$errors$R2,
pieColor = pieColor,
legend = TRUE, # <‑‑ leyenda automática
legend.cex = 0.8, # tamaño texto de la leyenda
layoutScale = c(0.8, 0.8), # encoge ligeramente la red
rescale = TRUE
)Figure 9. Networks of second Participants
first’s network
# --- first GROUP ---
# Convert data to matrix
TP_first_matrix <- as.matrix(TP_first)
# Estimate MGM model (defines number of variables, order of interactions, and input data)
invisible(capture.output({
TP_mgm_first <- mgm(
data = TP_first_matrix,
type = rep("g", length(TP_first)),
level = rep(1, length(TP_first)),
k = 2,
warnings = TRUE,
verbatim = TRUE
)
}))
# Predict explained variance (R²) and other node-level metrics
TP_predic_first <- predict(
TP_mgm_first,
TP_first,
error.continuous = "VarExpl"
)
# Plot the network with R² represented as pie charts on nodes
TP_glasso_predicted_first <- qgraph::qgraph(
cor_auto(TP_first),
graph = "glasso",
layout = TP_glasso_predicted$layout,
groups = group_list,
#color = group_colors,
palette = "colorblind",
nodeNames = item_labels,
labels = TRUE,
threshold = TRUE,
sampleSize = nrow(network_variables),
vsize = 5,
esize = 5,
asize = 6,
border.width = 1,
border.color = "black",
unCol = "black",
theme = "colorblind",
pie = SV_predic$errors$R2,
pieColor = pieColor,
legend = TRUE, # <‑‑ leyenda automática
legend.cex = 0.8, # tamaño texto de la leyenda
layoutScale = c(0.8, 0.8), # encoge ligeramente la red
rescale = TRUE
)Figure 10. Networks of first Participants
To evaluate the invariance of the centrality measures network across sexs.
# Compare centralities between networks
network_comparison <- compareCentrality(
TP_glasso_predicted_second,
TP_glasso_predicted_first,
include = "all",
legendName = "Centrality Measures by sex",
net1Name = "second",
net2Name = "first"
)
plot(network_comparison)Figure 11. Centrality comparison between second and first networks
set.seed(1234)
invisible(capture.output({
nct_test_centrality <- NCT(
EN_second,
EN_first,
it = 1000,
test.centrality = TRUE,
p.adjust.methods = "BH",
centrality = c("closeness", "betweenness", "strength", "expectedInfluence"),
)
}))
# Extract p-values and convert to data frame
centrality_pvals <- as.data.frame(nct_test_centrality$diffcen.pval)
# Round to 3 decimals (and adjust "< 0.001" if you want a more elegant output)
centrality_pvals_rounded <- round(centrality_pvals, 3)
# Optional: add row names as a column
centrality_pvals_rounded$Node <- rownames(centrality_pvals_rounded)
centrality_pvals_rounded <- centrality_pvals_rounded[, c("Node", setdiff(names(centrality_pvals_rounded), "Node"))]
row.names(centrality_pvals_rounded) <- 1:15
# Display as table
kable(centrality_pvals_rounded, caption = "Table 11. P-values for Centrality Differences (NCT Test)") %>%
kable_styling(full_width = FALSE, position = "center")| Node | closeness | betweenness | strength | expectedInfluence |
|---|---|---|---|---|
| age | 0.455 | 1.000 | 0.325 | 1.000 |
| AUTO | 0.325 | 1.000 | 0.312 | 1.000 |
| HETERO | 0.325 | 1.000 | 1.000 | 1.000 |
| TP | 0.455 | 0.356 | 0.280 | 0.738 |
| SP | 0.455 | 0.710 | 1.000 | 1.000 |
| RS | 0.455 | 1.000 | 1.000 | 0.312 |
| CIUS | 0.280 | 0.455 | 0.280 | 1.000 |
| IDS | 0.321 | 0.781 | 0.455 | 0.587 |
| AB | 0.587 | 1.000 | 1.000 | 1.000 |
| VB | 0.671 | 1.000 | 0.455 | 1.000 |
| VCB | 1.000 | 1.000 | 1.000 | 0.582 |
| ACB | 1.000 | 0.240 | 0.587 | 0.653 |
| AF | 0.280 | 0.587 | 0.653 | 0.280 |
| AA | 0.312 | 1.000 | 0.672 | 1.000 |
| APS | 0.312 | 1.000 | 0.487 | 0.710 |
The construction of the Bayesian network (BN) is carried out in two steps. Firstly, the estimation of the directed acyclic graph (DAG) is performed, and secondly, the BN model is fitted and validated with the study dataset.
For the DAG estimation phase, the PC Stable algorithm with no restrictions was employed. To ensure stability in the obtained DAG, a total of 200 bootstrap samples were drawn, and only the edges with a strength greater than 0.85 and a direction greater than 0.5 were retained in the final model.
set.seed(1234)
bootstr <- boot.strength(network_variables, R = 500, algorithm = "pc.stable")
avgnet <- averaged.network(bootstr, threshold = 0.85)
sp <- strength.plot(avgnet, bootstr, shape = "ellipse", render = FALSE)
TP_DAG_qgraph <- qgraph::qgraph(
sp,
layout = "spring",
groups = group_list,
#color = group_colors,
nodeNames = item_labels,
labels = TRUE,
sampleSize = nrow(network_variables),
vsize = 5,
esize = 5,
asize = 6,
border.width = 1,
border.color = "black",
unCol = "black",
theme = "colorblind",
legend = TRUE, # <‑‑ leyenda automática
legend.cex = 0.8, # tamaño texto de la leyenda
layoutScale = c(0.8, 0.8), # encoge ligeramente la red
rescale = TRUE
)Figure 12. Directed Acyclic Graph (DAG)
To fit and validate the BN model in the second phase of the process, the dataset was subdivided into 10 folds. A routine was implemented in which 90% of the folds were used to train the model, and the remaining 10% were used for testing. This cross-validation routine was repeated until all potential combinations were explore.
#Cross-validation (Check spatial cross-validation: blockCV)
set.seed(1234)
# Split data in 5 sets
kf <- dismo::kfold(nrow(network_variables), k = 10) # k-fold partitioning
models_fit <- c()
for (var in names(network_variables)) {
kfold_fit <- c()
aux_base <- data.frame('variable'=var)
if(is.numeric(network_variables[[var]])){
for(k in 1:10) {
test <- network_variables[kf == k, ]
train <- network_variables[kf != k, ]
training <- bn.fit(avgnet,train)
pred <- predict(training, node = var, data = test)
z <- data.frame( R2 = R2(pred, test[[var]]),
RMSE = RMSE(pred, test[[var]]),
MAE = MAE(pred, test[[var]]))
if(is.na(z$R2)){
z[is.na(z$R2),] <- 0
}
kfold_fit <- rbind(kfold_fit, z)
}
z <- apply(kfold_fit, 2, mean)
z <- cbind(aux_base,as.data.frame(t(z)))
models_fit <- rbind(models_fit, z)
}else{
for(k in 1:10) {
# Split data in test and train
test <- network_variables[kf == k, ]
train <- network_variables[kf != k, ]
training <- bn.fit(avgnet,train)
pred <- predict(training, node = var, data = test)
validation <- confusionMatrix(as.factor(pred),test[[var]])
z1 <- as.data.frame(t(validation$overall))
if(length(levels(network_variables[[var]])) > 2){
z2 <- as.data.frame(t(sapply(as.data.frame(validation$byClass),mean)))
}else{
z2 <- as.data.frame(t(sapply(validation$byClass,mean)))
}
z <- cbind(z1, z2)
kfold_fit <- rbind(kfold_fit, z)
}
z <- apply(kfold_fit, 2, mean)
z <- cbind(aux_base,z)
models_fit <- rbind(models_fit, z)
}
}
# Plot the network with R² represented as pie charts on nodes
TP_BN_qgraph <- qgraph::qgraph(
sp,
layout = "spring",
groups = group_list,
#color = group_colors,
palette = "colorblind",
nodeNames = item_labels,
labels = TRUE,
sampleSize = nrow(network_variables),
vsize = 5,
esize = 5,
asize = 6,
border.width = 1,
border.color = "black",
unCol = "black",
theme = "colorblind",
pie = models_fit$R2,
pieColor = pieColor,
legend = TRUE, # <‑‑ leyenda automática
legend.cex = 0.8, # tamaño texto de la leyenda
layoutScale = c(0.8, 0.8), # encoge ligeramente la red
rescale = TRUE
)Figure 13. Bayesian Network (BN) Model
# GGM - General
ggm_general <- SV_predic$errors %>%
transmute(
Variable = items,
R2_gen = round(R2, 3),
RMSE_gen = round(RMSE, 3)
)
# GGM - Female
ggm_female <- TP_predic_female$errors %>%
transmute(
R2_fem = round(R2, 3),
RMSE_fem = round(RMSE, 3)
)
# GGM - Male
ggm_male <- TP_predic_male$errors %>%
transmute(
R2_male = round(R2, 3),
RMSE_male = round(RMSE, 3)
)
# GGM - First
ggm_first <- TP_predic_first$errors %>%
transmute(
R2_first = round(R2, 3),
RMSE_first = round(RMSE, 3)
)
# GGM - Second
ggm_second <- TP_predic_second$errors %>%
transmute(
R2_second = round(R2, 3),
RMSE_second = round(RMSE, 3)
)
# BN model
bn <- models_fit %>%
transmute(
R2_bn = ifelse(R2 == 0, "-", round(R2, 3)),
RMSE_bn = ifelse(RMSE == 0, "-", round(RMSE, 3))
)
# Combine all into final table
final_table <- bind_cols(
ggm_general,
ggm_female,
ggm_male,
ggm_first,
ggm_second,
bn
)
# Set column names
names(final_table) <- c(
"Variable",
rep(c("R²", "RMSE"), times = 6) # General, Female, Male, First, Second, BN
)
# Generate final table
final_table %>%
kable(align = "lcccccccccccc", booktabs = TRUE,
caption = "Table 12. Explained variance and RMSE of partial-correlation network (GGM) and BN models") %>%
add_header_above(c(" " = 1,
"General" = 2,
"Female" = 2,
"Male" = 2,
"First Period" = 2,
"Second Period" = 2,
"DAG" = 2)) %>%
add_header_above(c(" " = 1,
"GGMs" = 10,
"BN" = 2)) %>%
kable_styling(full_width = FALSE, position = "center")| Variable | R² | RMSE | R² | RMSE | R² | RMSE | R² | RMSE | R² | RMSE | R² | RMSE |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| age | 0.166 | 0.912 | 0.180 | 0.905 | 0.146 | 0.923 | 0.102 | 0.946 | 0.057 | 0.970 | 0.116 | 1.178 |
| AUTO | 0.556 | 0.666 | 0.533 | 0.683 | 0.552 | 0.668 | 0.579 | 0.648 | 0.536 | 0.680 | 0.508 | 17.332 |
| HETERO | 0.483 | 0.718 | 0.433 | 0.752 | 0.507 | 0.701 | 0.473 | 0.725 | 0.488 | 0.715 |
|
|
| TP | 0.298 | 0.837 | 0.382 | 0.785 | 0.248 | 0.866 | 0.244 | 0.868 | 0.385 | 0.783 |
|
|
| SP | 0.452 | 0.740 | 0.436 | 0.750 | 0.461 | 0.733 | 0.371 | 0.792 | 0.460 | 0.734 | 0.415 | 0.617 |
| RS | 0.392 | 0.779 | 0.351 | 0.804 | 0.431 | 0.753 | 0.312 | 0.829 | 0.442 | 0.746 |
|
|
| CIUS | 0.431 | 0.754 | 0.536 | 0.681 | 0.335 | 0.814 | 0.458 | 0.736 | 0.434 | 0.751 | 0.414 | 0.525 |
| IDS | 0.373 | 0.791 | 0.414 | 0.765 | 0.323 | 0.822 | 0.412 | 0.766 | 0.336 | 0.814 | 0.18 | 0.748 |
| AB | 0.555 | 0.666 | 0.563 | 0.660 | 0.554 | 0.667 | 0.508 | 0.701 | 0.619 | 0.617 |
|
|
| VB | 0.528 | 0.687 | 0.601 | 0.631 | 0.507 | 0.701 | 0.509 | 0.700 | 0.561 | 0.661 | 0.476 | 0.652 |
| VCB | 0.511 | 0.699 | 0.412 | 0.766 | 0.578 | 0.648 | 0.438 | 0.749 | 0.586 | 0.643 |
|
|
| ACB | 0.573 | 0.653 | 0.451 | 0.740 | 0.637 | 0.602 | 0.465 | 0.730 | 0.666 | 0.577 | 0.518 | 0.406 |
| AF | 0.410 | 0.768 | 0.353 | 0.804 | 0.516 | 0.695 | 0.410 | 0.767 | 0.454 | 0.738 | 0.367 | 0.606 |
| AA | 0.447 | 0.743 | 0.512 | 0.697 | 0.390 | 0.780 | 0.463 | 0.732 | 0.408 | 0.768 | 0.431 | 0.607 |
| APS | 0.527 | 0.688 | 0.545 | 0.674 | 0.527 | 0.687 | 0.537 | 0.679 | 0.544 | 0.674 |
|
|